# libraries
library(tidyverse)
library(broom)
library(lme4)
library(here)
library(conflicted)
library(ggeffects)
library(Amelia)
library(broom.mixed)
library(ordinal)
library(patchwork)

# set working directory.
#setwd("")
print(getwd())

# resolve library conflicts
conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")


# read data
df = read_rds("r/replication/data/observational.rds") %>% 
  drop_na(vigilante_assault) %>% 
  mutate(vigilante_assault_f = factor(vigilante_assault, 
                                      levels = c(1, 2, 3), 
                                      labels = c("Neither approve nor understand", 
                                                 "Understand but not approve", 
                                                 "Approve")))


# clean dictionary
clean_dict = tribble(~dirty, ~clean, 
                     "severe_dum", "Victim of violent crime",
                     "police_patrols", "Frequency of police patrols", 
                     "help", "Police help citizens with problems", 
                     "assist_victims", "Police help victims", 
                     'police_latent', 'Attitudes towards police index', 
                     'broken_windows', 'Broken windows index',
                     'crime_problem', 'Crime is major community problem', 
                     'mano_dura', "The best way to reduce crime is with repression and mano dura")



# mixed effects ordinal
ordinal_mixed = clmm(vigilante_assault_f ~ severe_dum + age + woman + education + 
                       police_latent +  basic_expenses + indigenous_lang + 
                       km_to_cabecera + asset_count + war_victim + 
                       broken_windows_neigh + crime_problem + 
                       crime_better + (1 | municipio) +
                       (1 | departamento) + (1 | community), 
                      data = df)





# Mano dura: Figure 2, Table A.9 ---------------------------------------------------------------



# linear model
lm_dura = lmer(vigilante_assault ~ mano_dura + severe_dum +age + 
                 woman + education + 
                     police_latent +  basic_expenses + indigenous_lang + 
                     km_to_cabecera + asset_count + war_victim + 
                     broken_windows_neigh + crime_problem + 
                     crime_better + (1 | municipio) +
                     (1 | departamento) + (1 | community), 
                   data = df)



ordinal_dura = clmm(vigilante_assault_f ~ mano_dura + severe_dum + 
                      age + woman + education + 
                      police_latent +  basic_expenses + indigenous_lang + 
                       km_to_cabecera + asset_count + war_victim + 
                       broken_windows_neigh + crime_problem + 
                       crime_better + (1 | municipio) +
                       (1 | departamento) + (1 | community), 
                     data = df)





# marginal effects
mano_pred = ggpredict(ordinal_dura, terms = "mano_dura") %>% 
  as_tibble() %>% 
  mutate(var = "Mano dura") %>% 
  mutate(response.level = case_when(response.level == 3 ~ "Approve", 
                                    response.level == 2 ~ "Disapprove",
                                    response.level == 1 ~ "Understand but not approve")) %>% 
  mutate(label = str_replace_all(response.level, 
                                 pattern = "\\.", 
                                 replacement = " "))

police_pred = ggpredict(ordinal_mixed, terms="police_latent") %>% 
  as_tibble() %>% 
  mutate(var = "Police capacity index") %>% 
  mutate(response.level = case_when(response.level == 3 ~ "Approve", 
                                    response.level == 2 ~ "Disapprove",
                                    response.level == 1 ~ "Understand but not approve")) %>% 
  mutate(label = str_replace_all(response.level, 
                                 pattern = "\\.", 
                                 replacement = " "))


p1 = mano_pred %>% 
  filter(label == "Approve") %>% 
  mutate(var = "Support for harsh, authoritarian crackdown on crime") %>% 
  ggplot(aes(x = x, y = predicted, 
             ymin = conf.low, 
             ymax = conf.high)) + 
  geom_pointrange(fatten = 1, size = 3, fill = "white", shape = 21, 
                  color = "Coral3", alpha = .9) + 
  theme_light(base_family = "Fira Sans", 
              base_size = 14) + 
  theme(strip.background = element_rect(fill = "black"), 
        strip.text.y = element_text(size = 10), 
        strip.text = element_text(face = "bold"), 
        panel.grid.major = element_blank()) + 
  labs(x = "Strength of support for mano dura", y = "Pr(support vigilante killing)") + 
  scale_y_continuous(labels = scales::percent)

p2 = police_pred %>% 
  filter(label == "Approve") %>% 
  ggplot(aes(x = x, y = predicted, 
             ymin = conf.low, 
             ymax = conf.high)) + 
  geom_pointrange(fatten = 1, size = 3, fill = "white", shape = 21, 
                  color = "Coral3", alpha = .9) + 
  theme_light(base_family = "Fira Sans", 
              base_size = 14) + 
  theme(strip.background = element_rect(fill = "black"), 
        strip.text.y = element_text(size = 10), 
        strip.text = element_text(face = "bold"), 
        panel.grid.major = element_blank()) + 
  labs(x = "Local police capacity index", y = "Pr(support vigilante killing)") + 
  scale_y_continuous(labels = scales::percent) + 
  scale_x_continuous(breaks = c(-8, -1, 6), 
                     labels = c("Low", "Medium", "High"))


p1 / p2


ggsave(here("r", "replication", "figures", "mano-dura-obs.pdf"), device = cairo_pdf, width = 6, height = 8)



# need to make tables
cm <- c( '(Intercept)' = 'Constant', 
         'mano_dura' = "Mano Dura", 
         'severe_dum' = "Crime Victim",
         "police_latent" = "Police capacity index",
         "trust_pnc_f" = "Trust PNC",
         'age' = 'Age', 
         'war_victim' = "Civil war victim",
         'asset_count' = "Wealth (assets)", 
         'basic_expenses' = "Meets basic expenses", 
         "broken_windows_neigh" = "Broken windows index", 
         "crime_better" = "Safety has improved", 
         "crime_problem" = "Crime is major problem", 
         "education" = "Education", 
         "indigenous_lang" = "Primary indigenous lang.", 
         "km_to_cabecera" = "km to urban center", 
         "Neither approve nor understand|Understand but not approve" = "Disapprove --> Understand", 
         "Understand but not approve|Approve" = "Understand --> Approve")



modelsummary::msummary(list("Linear HLM" = lm_dura, 
              "Ordinal HLM" = ordinal_dura), 
         fmt = '%.3f', coef_map = cm, title = "Model output.", 
         label = "tab:obs", stars = FALSE,
         gof_omit = 'IC|Log|Adj|edf|REMLcrit', here("r", "replication", "figures", 
                                                    "tab-obs.tex"))
